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Abstract. - We theoretically study the excitation spectrum of confined macroscopic optical lattices in the 
Mott-insulating limit. For large systems, a fast numerical method is proposed to calculate the ground state filling and 
excitation energies. We introduce many-particle on-site energies capturing multi-band effects and discuss tunnelhng on 
a perturbative level using an effectively restricted Hilbert space. Results for small one-dimensional lattices obtained by 
this method are in good agreement with the exact multi-band diagonalization of the Hamiltonian. Spectral properties 
associated with the formation of regions with constant filling, so-called Mott shells, are investigated and interfaces 
between the shells with strong particle fluctuations are characterized by gapless local excitations. 



Introduction 

The equivalence of lattice sites is the foundation of solid state 
physics as it causes bands and lattice-periodic wave functions. 
Interesting physics emerges when the band is bend caused by 
the inhomogeneity of the system, e.g. in clusters, where bulk 
and surface electrons show a different physical behaviour. In 
optical lattices, the inhomogeneity comes naturally due to the 
finite waist of the lattice-establishing laser beams or due to an 
additional dipole trap. It has been shown in a series of pio- 
neering work [1-3] that the bosonic Mott-insulator phase can 
be realized in optical lattices. Further measurements [4-6] 
have demonstrated that the situation is more subtle and Mott- 
insulator shells appear, i.e. plateaus with constant filling fac- 
tors descending in integer steps from the centre of the trap. 
The existence of superfluid regions between the Mott shells has 
been theoretically predicted [7-12]. The system has previously 
been studied numerically, using quantum Monte Carlo [7-9] 
and DMRG [10], and analytically in Refs. [11, 12]. In addition 
to condensed matter aspects, the coexistence of compressible 
and incompressible regions has important implications on adi- 
abatic heating in optical lattices [13]. 

We present a multi-band exact diagonalization study of small 
systems exploring the exact excitation spectrum and the pre- 
cursor of shell formation. Based on this calculation we use a 
numerical method, which in the limit of deep lattices allows 
us to obtain an approximated energy spectrum and the exact 
occupation numbers. Thereby, we incorporate orbital changes 
by introducing a particle-number-dependent on-site interaction. 
The formation of shells, local excitation gaps and particle fluc- 
tuations are discussed reflecting the strong inhomogeneity of 



the system. Our approach is suitable for optical lattices with 
millions of atoms in arbitrary spatial dimensions and allows a 
perturbative treatment of tunnelling. Compared to algorithms 
such as quantum Monte Carlo [7-9] and DMRG [10], this tech- 
nique gives numerically inexpensive results for large 3D lattice 
systems in a specific parameter regime. We organize this paper 
as follows. Subsequent to the introduction of the system, we 
present results obtained by means of exact diagonalization fol- 
lowed by the numerical approach for vanishing tunnelling and 
a comparison of both results. Finally, we discuss consequences 
for two- and three-dimensional optical lattices including non- 
zero tunnelling. 

The interaction of the ultracold bosonic atoms with mass m 
is modelled by a contact potential ^^(r — r') with g = ^^«s 
and the s-wave scattering length a^. Using the bosonic field 
operator ^(r), the Hamiltonian including the repulsive two- 
particle interaction reads 



2m 



+ Fp(r) + yc(r) 



^(r) 



(1) 



where the periodic potential Vp of the optical lattice is given by 
Vb,cc cos^(7rx/a) + Vb,2y cos^ (ttt// a ) H-Vb,^ cos^(7rz/a) with the 
lattice constant a. The atoms experience an additional confine- 
ment potential Vc caused by the finite Gaussian beam waist Wq 
of the laser beams and/or an additional dipole trap with the fre- 
quency ujd- Using a harmonic approach, the confining potential 
in x-direction is given by 



Vc = Vh 



,2' 



(2) 
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Fig. 1: (a) The periodic potential Vp is superposed by a harmonic potential Vc leading to site offsets e^. The single-particle spectrum for 
Vh — 0.06^H containing discrete bands is plotted within the potential, (b) The many-particle energy spectrum (91 states) obtained by exact 
diagonalization (upper half) and for vanishing tunnelling (lower half) as a function of the harmonic confinement strength Vh. Additional states 
in the J ^ spectrum are plotted in grey. The energy scale is relative to the ground state energy and the gap to the first excited state is 
depicted in grey. 



where Vh = Im^^a^ with u^l^ « + 
Exact diagonalization 

For the exact diagonalization, the potential is truncated to a ID 
lattice with rux = 10 and rriy = = 1 sites by adding 
a smooth boundary [14]. Exemplarily, we have chosen ^^Rb 
atoms with = lOOao, where ao is the Bohr radius, = 10 
particles, a lattice constant a = 515 nm, and a transversal 
lattice depth Vo,y = Vo,z = ^OEr, where Er = is 
the recoil energy. The lattice depth in the longitudinal direc- 
tion is Vq^x = '^OEr, leading to a system deep within the 
Mott-insulator phase for a vanishing harmonic confinement 
(Vh = 0). The potential of the ID lattice is shown in figure la, 
depicting the energy offsets of the modeled chain, which are 
labeled eo = for the two equivalent central sites and 64 for the 
outermost sites. The parabolic shape of the confinement leads 
to large energy offsets of outer sites, although the offsets are 
small for central ones. At = 0, the single-particle spectrum 
shows three discrete bands that comprise ten delocalized Bloch 
states each (corresponding to the number of sites), where the 
width of the lowest band is 4 J ?^ 2.4 x 10~^ Er [15]. For larger 
values of Vh the 'band' width is given basically by the offset of 
the outermost sites (see figure la). At > OAEr, the outer- 
most site offset 64 > SEr causes the lowest two single-particle 
bands to overlap. 

To calculate the many-particle spectrum by means of exact 
diagonalization, we limit the many-particle basis to 500 000 
Fock states with lowest energy and freeze out the orbital de- 
grees of freedom in y- and z-direction. Afterwards, we calcu- 
late the matrix elements of (1) and obtain the 91 lowest eigen- 
values (the ground state and the 90 states of the first band for 
Vh 0). The exact diagonalization method includes parti- 
cle correlations and the admixture of higher bands, i.e., orbital 
changes taking place for higher fillings [14], but is strongly lim- 
ited to systems with few lattice sites. The many-particle spec- 
trum is plotted in figure lb (upper half) relative to the ground 
state energy £^0 as a function of the harmonic confinement Vh 
using a logarithmic scale. For convenience, the band gap be- 
tween the ground state and the first excited state is depicted in 



grey. 

For vanishing confinement (Vh = 0), the spectrum shows 
that the first excited band is gapped from the ground state by the 
on-site interaction energy of two particles U ^ 0.6Er, which 
corresponds to the gap of the macroscopic Mott-insulator 
phase. The atoms are strongly localized on single lattice sites 
[14], so that all lattice sites are occupied by exactly one particle 
per site. Increasing the harmonic confinement leads to abrupt 
crossovers to states with higher integer occupation numbers, 
i.e. finite size correspondents to Mott shell configurations. In 
the many-particle spectrum these crossover points are reflected 
by vanishing energy gaps. In between these points lobe-like 
energy gaps can be observed, where the lobes correspond to 
the occupation number configurations 1-1-1-2-2-1-1-1, 1-2-2- 
2-2-1, 2-3-3-2, 1-4-4-1, and finally for Vh > IEr a states with 
5 particles on the two central sites [16]. For Vh < O.OSEr, 
the lattice is commensurately filled with one particle per site. 
In this region, the excited band broadens with increasing Vh, 
leading to a slowly decreasing energy gap. At Vh ~ 0.03Er, 
where the energy gap vanishes, a multitude of low-lying exci- 
tations is possible leading to a compressible system. At this 
point, the double occupation of the central sites, which corre- 
sponds to the on-site energy U, becomes energetically prefer- 
able to the occupation of the outermost sites with 64 = 20Vh. 
Analogous situations exist at the other crossovers, so that the 
shell structure is dominated by the ratio of on-site interaction 
ni{ni-i) jj individual site offsets e^. The critical behaviour 
separating different filling configurations manifests itself in the 
disappearance of the excitation gap. 

Classical approach and orbital degrees of freedom 

The nature of the excited states, however, is more complicated 
and the energy spectrum, which contains important physics, is 
rather complex. In the following, we therefore consider the 
'classical' case of vanishing tunnelling ( J ^ 0) which is valid 
deep within the Mott-insulator phase. Using this approach, we 
show that the basic features of the spectrum can be uncovered. 
Furthermore, larger systems can be studied, for which finite 
tunnelling can be reintroduced perturbatively in a second step. 
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with this method but their positions are shifted. This effect 
is additionally enhanced for higher filling factors. To study 
this influence more quantitatively, we calculate the on-site few- 
particle interaction with exact diagonalization. Since the lattice 
depth is sufficiently deep, it is valid to approximate the lowest- 
band Wannier function W(^(v) by the wave function of a single 
sinusoidal lattice site using hard boundary conditions. The cor- 
relationless Hubbard on-site interaction U = g J d^r \ wo{r)\^ 
does not incorporate that particles tend to avoid each other for 
repulsive interaction. This causes broadening of the particle 
density and was addressed experimentally in [5]. Using the 
correct many-particle wave function for n particles, the ex- 
pectation value of the on-site interaction becomes 



Fig. 2: A comparison of Hubbard on-site interaction U, the many- 
particle interaction energy u{n), and the total on-site energy e(n) per 
atom pair for Vb,i = Vb and i — x, y, z. The single-site many- 
particle wave function is determined by exact diagonalization for n = 
2 particles (highest energy) to n = 5 particles (lowest energy). 

dealing with a drastically reduced basis set. For J 0, the 
truncated Bose-Hubbard Hamiltonian is given by 

H^Y.'^'^^U^e.n. (3) 

i 

and the localized occupation number basis |ni, n2, um) is 
an eigenbasis of the Hamiltonian, where M = mxruymz de- 
notes the number of sites. In principle, finding the ground state 
requires to calculate the total energy E = (H) of all possible 
basis states. The efficiency of this method is very limited due 
to the huge number of basis states ^n\{m-i)\ large lattices, 
where N is the number of particles. Following [11, 12], in the 
local density approximation an effective local chemical poten- 
tial /i(r) = /i — Vc{r) can be introduced, and the ground state 
is then constructed via filling each lattice site up to the local 
chemical potential (1 separately. Using the continuous limit, 
the chemical potential /i(A^) is calculated analytically in [11] 
and numerically in [12]. The continuous limit is, however, 
only applicable to smooth confining potentials. In general, the 
self-consistent determination of /i(A/') is tedious. We therefore 
suggest the following iterative algorithm for the solution of (3) 
directly in the microcanonical ensemble with a fixed total par- 
ticle number N = n^. It allows us also to construct the 
lowest excited states of the system and therefore a comparison 
with figure lb (upper half). Starting with an empty lattice, the 
N particles are added successively to that site j of the lattice, 
where the expense of energy fi^ = UjU -\- Cj is presently mini- 
mal. This procedure gives the lowest energy occupation for any 
particle number A^, since /i^ > and all sites are uncorrected. 
The complexity of this algorithm is given by 0{NM), since 
each step adds one particle to one of M possible sites. This 
allows the calculation of the exact occupation numbers for a 
million particles, e.g., N = M = 100^, within seconds on an 
ordinary desktop computer. This method is therefore consider- 
able useful for the design and interpretation of experiments. 

Compared with the exact diagonalization results for the ID 
lattice, the crossover points can, in principle, be reproduced 



u{n) = -^-^(V,t(r)^t(r)^(r)V;(r)), (4) 

which is normalized to the interaction energy per atom pair. 
The results that are depicted in figure 2 show that u{n) de- 
termined by exact diagonalization deviates strongly from the 
Hubbard on-site interaction U. As expected, the interaction en- 
ergy u{n) decreases with an increasing number of particles per 
site. However, not only u{n) changes when the modification 
of wavefunctions is taken into account. In fact, the admixture 
of correlated states and the broadening of the density change 
the expectation value of on-site kinetic and potential energy. 
The total on-site energy is the eigenvalue of the many-particle 
Schrodinger equation Hi = restricted to single-site 

wave functions and using the full Hamiltonian (1) with Vh = 0. 
The normalized on-site energies e(n) = 2(£'^ — E^)/n{n — 1) 
are plotted in figure 2, where E^ is the energy of the non- 
interacting system. It shows, that the total interaction lies en- 
ergetically between U and u{n). Nevertheless, the deviation 
of e(2) from the Hubbard U is large. Using U = e(2), the 
crossovers are still shifted noticeably for higher fillings (large 
Vh) comparing with the results for the ID lattice (upper half 
of figure lb). Therefore, we incorporate in our calculations the 
correct total interaction energy e{nj) for a single site with the 
occupation rij . Please note that the optimization problem re- 
mains the same if substituting the energy /i^ = '^jU + ej for 
adding one particle by /i^ = e{nj-\-l)—e{nj)-\-ej. In the lower 
half of figure lb the energy spectrum for vanishing tunnelling 
(J 0) is shown using the corrected values of the on-site in- 
teraction. In this case, the crossover energies corresponding to 
the vanishing gaps for J ^ are in good agreement with the 
exact diagonalization of the ID lattice. This result shows, in 
general, that the introduced filling-dependent on-site interac- 
tion e(n) is appropriate to describe effects arising from orbital 
changes. The small remaining shift of the crossover energies in 
figure lb is due to the 'classical' treatment of the states in our 
approach. 

The energy gap can be obtained by removing one particle 
from site j and adding it to site k ^ j. Finding the minimum 
excitation energy AEj^k for all possible j and k has in gen- 
eral the complexity O(M^). The excitation energy is given by 

AEj^k = y^j^ + Mfc ' where 

/i+ = e{nj + 1) - e{nj) + e^- > (5) 
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for adding a particle on site j and 

11^ = e{nk - 1) - e(n/e) - < (6) 

for removing a particle from site k. Thus, it is sufficient to 
minimize ji^ and /i^ separately, which reduces the complex- 
ity to linear order in M. Finding the next excited state is more 
complicated, since this state may be an excitation of the ground 
state but also of the first excited state. For the calculation of the 
spectrum, a slightly modified Dijkstra algorithm can be used, 
where one considers the ground state as the node of a graph. 
This node is expanded according to all possible excitations. It- 
eratively, the node with the minimum energy, which is not ex- 
panded yet, is expanded. Consequently, the list of expanded 
nodes represents the states with the lowest energy. It is neces- 
sary, however, to check whether a created excited state is al- 
ready a node in the graph, which is a major contribution to the 
complexity of the Dijkstra algorithm, but can numerically be 
highly optimized. Using this procedure, the resulting energy 
spectrum (lower half of figure lb) reproduces well the basic 
features, i.e. the band gap, the overall shape and the density of 
states, of the exact calculation (upper half). Because of inter- 
actions, the degeneracy of states is lifted in the exact spectrum, 
so that the actual energies of many states are shifted. 

Macroscopic lattices 

Transferring the above results to macroscopic lattices is not 
straightforward. Enlarging the system's size causes the differ- 
ences in offset energies of neighbouring sites Cj — Ck to de- 
crease substantially when keeping the offset of the outermost 
sites fixed. This causes the width and the height of the energy 
lobes in the spectrum to decrease because a huge number of 
configurations become possible when increasing the number of 
particles and lattice sites. This process is drastically enhanced 
in 2D and 3D lattices, where practically the band gap vanishes 
for all confinement strengths Vh- The only exception is the real 
Mott-insulator phase, where the outermost sites have an offset 
smaller than U. In fact, the excitation spectrum of the total 
system becomes more or less continuous. This might appear 
contradictory to the insulating property at first glance but the 
system is inhomogeneous and only some of the atoms can per- 
form gapless excitations. Therefore, local properties are more 
suited to describe the system, such as particle number fluctu- 
ations and the local excitation spectrum. It has already been 
pointed out [8-1 1] that the regions with fixed occupation num- 
ber, the Mott-insulator shells, are surrounded by compressible 
shells for non-negligible tunnelling. 

Exemplarily, the occupation number distribution ni for TV = 
10^ particles in a 3D lattice is depicted in figure 3a. It shows the 
expected shell structure with filling factors five (in the centre) to 
one. An advantage of the presented numerical algorithm is that 
it is also applicable to rapidly varying confinement potentials. 
It can thus be used to tailor more sophisticated shell configura- 
tions for experiments. Since the on-site interaction energy U is 
relatively small compared to the depth of the lattice wells, the 
local filling factors can be adjusted without a stronger pertur- 
bation of the 3D lattice. Additional laser beams or magnetic 
fields can thus be used to obtain complex spatial distributions 
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Fig. 3: (a) Cut through the occupation number distribution m of a 3D 
lattice {N = 10^ J = 0, 14 = 7.5 x 10"^^^ and Vo = 2^Er) 
showing fillings = 1 (outer shell) to = 5 (inner shell), (b) 
Particle fluctuations An^ due to finite tunnelling {J — 10~^Er ^ 
2 X 10~^U) appear at the surfaces of the shells {z — plane). 



of atoms with specific filling factors. Figure 4a shows the oc- 
cupation number distribution for an added periodic potential in 
x-direction Vc = VhX^ / o? co^^{itx/ a') with the periodic- 
ity a' = 3a motivated by [17]. The superlattice structure leads 
to alternating local fillings. As a second example, an added 
hat-shaped (— |x|) potential causes two separate atomic clouds, 
shown in figure 4b, that are in touch with each other at the ori- 
gin. 

Single particle excitations and finite tunnelling 

At the outer surface of each shell, almost gapless excitations are 
possible via the hopping of a particle to another site, whereas 
the inner surface can easily absorb particles. This is shown 
quantitatively in figure 5a (upper half), where A£^~ = + 
minj(/i^) is the minimal energy for removing one particle 
from site i and adding it to another site j. The minimal en- 
ergy for adding one particle to site i (and removing it from 
site j) is denoted as /S.Ef = /j,^ -\- minj(/i~) (lower half of 
figure 5 a). The local excitation gap vanishes at the outer and 
inner surfaces and increases strongly within the shells. The 
excitation energies AEf are explicitly shown for sites with 
y = z = in figure 5b. It reflects basically the harmonic 
shape of the confinement, so that /^E^ and A£^~ are posi- 
tive and negative parables, respectively, subtracted by integer 
multiples of the on-site interaction. The situation changes dras- 
tically, if accounting only for excitations to nearest neighbours 
as shown in figures 5c and 5d. Low energy excitations (either 
A£;+ or AE-) can occur exclusively on sites directly at the 
boundary between different shells. Within the bulk the single- 
particle excitation gap is wide and nearly constant. Therefore, 
the nearest-neighbour tunnelling is, in general, strongly sup- 
pressed on these sites, where the system is a good insulator. 

In the following, we use a perturbative approach to obtain 
particle fluctuations for finite tunnelling J. For each site i the 
subsystem containing the site i and all neighbouring sites is 
considered. The diagonalization of this subsystem with finite 
tunnelling — j) ^H^ws an approximative calculation 

of the particle fluctuations An^ = \/ (^f ) — (^i)^. For the ex- 
ample above, weak tunnelling causes finite fluctuations An^ at 
the boundaries of the shells, which is shown in figures 3b and 
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Fig. 4: (a) An additional standing wave with periodicity 3a and am- 
plitude U causes alternating fillings in x direction (N = 3.5 x 10^). 
(b) An added — 5014|x| potential leads to two separate atom spheres 
{N = 1.5 X 10^). The fluctuations Am for J = 10"^^^ and z = 
are shown at the bottom (see figure 3 for the colour map and other 
parameters). 




Fig. 5: The minimal excitation energy in the z = plane (a) and for 
2/ = z = (b) for particles hopping from site i (E~) and to site i 
(E^). The excitation energies only accounting for nearest-neighbour 
hopping in the z = plane (c) and for ^ = 2; = (d) including 
particle fluctuations An^ for J = 5 x 10~^ Er. 

4 for J = and in figure 5d for J = 5 x IO'^Er. 

Please note, that due to the finite cell size of the lattice, the re- 
sults are not completely spherically synmietric. In accordance 
with the nearest-neighbour excitation energies, the fluctuations 
affect only sites directly at the surfaces. Because of the larger 
slope of the confinement and the lower filling, the tunnelling 
decreases for outer shells. Within the bulk of the Mott shells 
small particle fluctuations can be observed due to the finite tun- 
nelling. The fluctuations on the surface however enhance the 
tunnelling of atoms next to the surface, which is not covered by 
this perturbative approach and would require a self-consistent 
calculation. The narrow energy gaps of particles close to the 
surface can, in principle, cause near-resonant tunnelling to non- 
nearest neighbours, so-called variable range hopping [18]. Due 
to the relatively strong Vh (compared with J) and the regularity 
of the potential, only sites on the surface, which are dominated 
by nearest-neighbour hopping, have suitably low excitation en- 
ergies. 

The presented method provides the lowest energetic states 
(including particle-hole excitations) and is suitable to restrict 
the Hilbert space effectively. The obtained states can be used 
as a starting point for exact diagonalization, quantum Monte- 



Carlo, and dynamical mean-field calculations including tun- 
nelling and finite temperatures. In particular, the diagonaliza- 
tion of the ID lattice (figure lb) could be transformed with 
much fewer basis states. 

Conclusions 

In summary, we have studied the excitation spectrum and ex- 
act site occupation numbers for confined optical lattices deep 
within the Mott insulator regime. In good agreement with ex- 
act diagonalization for small ID lattices, a numerical method 
was presented that allows for negligible tunnelling the exact 
treatment of macroscopic optical lattices with arbitrary shape 
of the confining potential. Adding slowly varying potentials 
to the optical lattice can give rise to complex filling structures. 
We have calculated the numerically exact many-particle on-site 
energies and have shown that introducing a filling factor de- 
pending on-site interaction can be incorporated to cover orbital 
changes. For small systems, the many-particle spectrum con- 
tains lobes, whereas for macroscopic systems nearly gapless 
excitations are always possible at the boundaries of the Mott 
shells. Within a given Mott shell the local excitation energy 
varies strongly leading to compressible sites close to the surface 
and incompressible inner sites. A perturbative treatment for fi- 
nite tunnelling shows strong particle fluctuations at the bound- 
aries between the shells, where the Mott-insulator gap vanishes. 
Finally, the presented method can serve as a reduction scheme 
for the Hilbert space for further numerical treatments. 
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